%-------------------------------------------------------------------------;
% Horizon Bias and the Term Structure of Equity Returns;
% This code replicates Table 5;
%-------------------------------------------------------------------------;

clear all;
clc;
close all;

%-----------------------------;
% Import term structure data;
%-----------------------------;

filename='Term_structure_variables';
data_ts=xlsread(filename);

b=size(data_ts,1); %Length of the data;

Div_ret=data_ts(2:end,6); %BBK extended;
SP_ret=data_ts(2:end,4);

Div_ret=log(1+Div_ret);
SP_ret=log(1+SP_ret);
                                                 
etp=SP_ret-Div_ret;

SP_price=data_ts(1:end,1);
SP_Div=data_ts(1:end,3);
SP_strip_18=data_ts(1:end,5);

SP_pd=log(SP_price./SP_Div);
SP_pd_strip=log(SP_strip_18./SP_Div);

SP_pd=SP_pd(1:b-12);
SP_pd_strip=SP_pd_strip(1:b-12);

%-------------------------------;
% Import explanatory variables;
%-------------------------------;

filename='Explanatory_variables';
data_ev=xlsread(filename);

data_ev=data_ev(1:b-12,:);

year=data_ev(1:end,1);
month=data_ev(1:end,2);

Ltg=data_ev(1:end,3);   
Stg=data_ev(1:end,4);   

%Take logs;
Ltg=log(1+Ltg);
Stg=log(1+Stg);

%Unaffiliated;
Ltg_un=data_ev(1:end,5);
Stg_un=data_ev(1:end,6); 

Ltg_un=log(1+Ltg_un);
Stg_un=log(1+Stg_un);

%Skewness;
skew1=data_ev(1:end,7);
skew1=(skew1-mean(skew1))./std(skew1);

skew5=data_ev(1:end,8);
skew5=(skew5-mean(skew5))./std(skew5);
Diff_skew=skew5-skew1;

%Uncertainty;
uncert1=data_ev(1:end,9);
uncert1=(uncert1-mean(uncert1))./std(uncert1);

uncert5=data_ev(1:end,10);
uncert5=(uncert5-mean(uncert5))./std(uncert5);
Diff_uncert=uncert5-uncert1;

%Other variables capturing business cycle;
rec=data_ev(1:end,11);
default=data_ev(1:end,12);
term=data_ev(1:end,13);
cay=data_ev(1:end,14);

%Horizon bias;
acorr=autocorr(Stg);
theta_0=mean(Stg);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg;

Stg_impl=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB=Ltg-Stg_impl;

%Horizon bias - unaffiliated;

acorr=autocorr(Stg_un);
theta_0=mean(Stg_un);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_un;

Stg_impl_un=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB_un=Ltg_un-Stg_impl_un;

%-------------------------;
% Predictive regressions;
%-------------------------;

for h=6:6:60

%Number of lags for overlapping observations;
lags=h-1;

%Sum returns over h periods;
ret_SP=filter(ones(h,1),1,SP_ret);
ret_Div=filter(ones(h,1),1,Div_ret);

ret_SP=ret_SP(h:end);
ret_Div=ret_Div(h:end);

Ltg_h=Ltg(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
Stg_h=Stg(1:min(size(Stg,1),(size(Stg,1)-h+12)));

Ltg_h_un=Ltg_un(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
Stg_h_un=Stg_un(1:min(size(Stg,1),(size(Stg,1)-h+12)));

ret_SP=ret_SP(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
ret_Div=ret_Div(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

SP_pd=SP_pd(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
SP_pd_strip=SP_pd_strip(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

HB=HB(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
HB_un=HB_un(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

rec=rec(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
default=default(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
term=term(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
cay=cay(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

skew1=skew1(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
skew5=skew5(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

%Dependent variable;
y=(ret_SP-ret_Div);

%HB + controls + error + NBER;
x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB SP_pd term default cay error rec];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Tab_5(1:14,h)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1);b(8,1);tb(8,1)];
Tab_5(15:16,h)=[length(error);AR2];
end;

%------------------------------------------------;
%Table 5;
%------------------------------------------------;

Table5=Tab_5(:,[6 12 18 36 60])
